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ABSTRACT 

A simple method is presented for the rapid simulation of statistically-isotropic non-Gaussian 
maps of CMB temperature fluctuations with a given power spectrum and analytically- 
calculable bispectrum and higher-order polyspectra. The nth-order correlators of the pixel 
values may also be calculated analytically. The cumulants of the simulated map may be used 
to obtain an expression for the probability density function of the pixel temperatures. The 
statistical properties of the simulated map are determined by the univariate non-Gaussian 
distribution from which pixel values are drawn independently in the first stage of the simu- 
lation process. We illustrate the method using a non-Gaussian distribution derived from the 
wavefunctions of the harmonic oscillator. The basic simulation method is easily extended to 
produce non-Gaussian maps with a given power spectrum and diagonal bispectrum. 
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1 INTRODUCTION 

The study of non-Gaussianity of Cosmic Microwave Background 
(CMB) fluctuations is of major importance in understanding the 
processes responsible for generating the fluctuations and in assess- 
ing the contribution of foreground astrophysical processes and in- 
strumental effects to observations of the CMB. 

Single-field inflationary scenarios predict, in gen- 
eral, that C MB fluctuations are very nearly Gaussian (see 
e.g. iBartolo et alj 120041 for a recent review), if one as- 
sumes that the sub-Hubble-sca le quantum fluctuations 
start off in the ground state ( Contaldi. Bean & Magueiio 
[ 19991: iMartin. Riazuelo & SakellarioudovJ ~ l2000r. 

iGaneui. Martin & Sakellari oudoul 120021) . In such models, non- 
Gaussianity produced during inflation arises predominantly from 
the non-linear nature of gravitational interactions rather than from 
self-interaction of the fluctuations of the inflaton field. Typically the 
level of non-Gaussianity is suppresse d by the first-or der slow-roll 
parameters (Acrjuaviva et al. 2003; Maldacena 2003). Subsequent 
non-linear processing of the primordial fluctuations to second- 
order in perturbation theory has been shown to amplify the tiny 
primordial non-Gaussianity to a level on the last-scattering surface 
that m ay be detectable with future CMB surveys JBartolo et alJ 
2004). Furthermore, second-order radiative transfer effects, 
such as gravitational lensing of the CMB, should produce a 
detec table level of non-Gaussianity in the CMB (Zaldarriaaa 
2000). Larger levels of non-Gaussianity can be produced in 
inflation models with multiple scalar fields. Examples include the 
curvaton (e.g.lLvt h & Wands 2002) and t he inho mogeneous reheat- 
ing (e.g. iDvali. Gruzinov & Zaldarriaeal 120041) scenarios. Finally 



models which include topological defects also produce signifi- 
cantly non-Gaussian fluctuations iGangui. Pogosian & Winitzka 
2002). It is also the case that inevitable contaminants, such as 
discrete radio sources, Galactic emission and systematic instru- 
mental effects leave non-Gaussian signatures on CMB maps. Thus 
non-Gaussianity tests are of fundamental importance both for 
probing inflation physics and for isolating systematic effects. 

In order to investigate one's ability to detect and recover non- 
Gaussian signals, it is useful to generate non-Gaussian maps with 
known statistical properties, to which putative analysis methods 
may be applied. In particular, in many applications, it is desir- 
able that the simulated non-Gaussian map is statistically isotropic, 
with a prescribed power spectrum (or 2-point correlation func- 
tion). The generation of such non-Gaussian CM B maps is, how- 
ever, a surprisingly difficult task (see, for example, |Vtoet_al. 2001, 
2002). Moreover, it is often useful for the non-Gaussian map 
also to have a known (or prescribed) bispectrum and one-point 
marginal probability density function. Several different techniques 
have been proposed to add ress various subsets of these re quire- 
ments l Contaldi & Magueiio 2001; Martinez-Gonzalez et al. 2002; 
Komatsu et al. 2003; Liguori et al. 2003), but each carries a con- 
siderable computational cost. The aim of this paper is to present 
a simple, fast technique for simulating statistically-isotropic non- 
Gaussian CMB maps with a prescribed power spectrum, for 
which one can calculate analytically the bispectrum, higher-order 
polyspectra, nth order pixel correlators and the one-dimensional 
marginalised distribution (in terms of its cumulants). Moreover, 
our basic simulation method is easily extended to produce non- 
Gaussian maps with a given power spectrum and diagonal bispec- 
trum. 
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The problem of simulating non-Gaussian CMB maps can be 
formalised as follows. A real, random scalar field T(x) can be 
defined as a collection of random variables, one at each point 
x = (x\,X2,.. . ,jc„) in the ^-dimensional space (clearly n — 2 for 
CMB maps on the sphere or flat-patches of sky). Thus, for each po- 
sition x, T(x) = t, where t is a scalar random variable with a one- 
dimensional (marginal) probability density function (PDF) fr(t). 
For a random field that is statistically homogeneous, frit) is the 
same at all points in the space. The main difficulty in the numerical 
simulation of a generic random field is that, in general, given two 
arbitrary positions x\ and X2, the quantities T(x\) and T(x 2 ) are 
not independent. In particular, it is often desirable for T(x) to have 
a prescribed 2-point covariance function 

% T {x u x 1 ) = {T{x x )T(x 2 )). 

In the case of a statistically homogeneous random field in Eu- 
clidean space, the covariance function depends only on r = x\ —x 2 - 
If the field is also isotropic then the dependence is only on x = |r|. 
As discussed bv lVio et all200ll l2002. most methods for sim- 
ulating non-Gaussian maps with a prescribed 2-point covariance 
function [and a prescribed marginal PDF fr(t)] are based on first 
generating a zero-mean, unit- variance Gaussian random field G(x), 
with an appropriate covariance structure ^cM- One then performs 
the mapping transformation G(x) — » T(x) according to 

T(x)=h[G(x)}, 

where h represents an appropriate function. The usefulness of this 
approach derives from the fact that there exist explicit analytical 
(but complicated) formulae for the marginal PDF, fr(t), and the 
covariance function, ^r(T), of the transformed field in terms of 
the covariance function, ^c(x), °f the original Gaussian field and 
the mapping function h. In particular, we note that the formula for 
^7-(x) takes the form of a double integral of a function depending 
on both %c(i) and h. There exist only a few mapping functions h 
for which fr{t) and ^r(t) ma y be calculated analytically. A still 
smaller subset of these cases allows the resulting expressions to be 
inverted analytically to obtain the required functions %c(i) an d h 
to be used in simulating the non-Gaussian map. In general, one has 
to resort to numerical methods to invert the general formulae for 
fr(t) and |r(x) and this can be computationally very costly. 

As mentioned above, it is often desirable for the simulated 
non-Gaussian map also to have a known (or prescribed) bispec- 
trum. The method outlined above has not been extended to this 
case, and any such generalisation is likely to be extremely com- 
putationally demanding. An alternative method for generating non- 
Gaussian maps with a prescribe d power spec t rum a nd bispectrum 
has been suggested by Contaldi & Maeueiio (2001), although, in 
general, the marginal distribution fj (t) of the resulting map cannot 
be obtained analytically. The method is based on choosing some 
one-dimensional non-Gaussian PDF, from which the real and imag- 
inary parts of (some subset of) the spherical harmonic coefficients 
a{ m of the map are drawn independently (the remaining coefficients 
being drawn from a Gaussian PDF). However, statistical isotropy 
imposes 'selection rules' upon correlators of the ci( m coefficients. 
For example, the third-order correlators must satisfy 
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where (■ ■ •) denotes the Wigner 3/ symbol and B^^ are the bis- 
pectrum coefficients. Hence the map corresponding to the drawn 
ag m values is not only non-Gaussian but also anisotropic since 
all its third-order correlators are zero, except for a subset of the 



form {cq ) for those ci( m drawn from the non-Gaussian PDF. It is 
therefore necessary first to produce an ensemble of non-Gaussian, 
anisotropic maps and then create an isotropic ensemble by applying 
a random rotation to each realisation. Contaldi & Masueiio 1 2001 ) 
show that these random rotations produce the necessary correla- 
tions between the ag m coefficients to ensure isotropy of the ensem- 
ble, but the method clearly requires significant computation. 

In this paper, we discuss a very simple and computationally 
fast method for simulating non-Gaussian maps that are, by con- 
struction, statistically isotropic, and for which numerous statistical 
properties may be calculated analytically. In particular, it is possi- 
ble to produce a map with a prescribed power spectrum for which 
one can obtain simple analytical expressions for the bispectrum and 
higher-order polyspectra, and nth-order correlators of the pixel val- 
ues. One may also calculate the cumulants of the map, which may 
be used to obtain the one-dimensional marginalised distribution. A 
simple extension of the method allows for the simulation of maps 
with a prescribed power spectrum and diagonal bispectrum. In Sec- 
tion|2| we discuss our method for simulating non-Gaussian maps, 
the statistical properties of which are presented in Section In 
Section[4]we extend the method to allow simulation of maps with 
prescribed power spectrum and diagonal bispectrum. Finally, our 
conclusions are presented in Section|5| 



2 SIMULATION METHOD 

Since our goal is the simulation of CMB maps for use in later anal- 
ysis, it is convenient at the outset to divide the celestial sphere into 
pixels labelled by p = 1,2,... , /V p j x . For simplicity, we also assume 
an equal-area pixelisation, so that each pixel subtends the same 
solid angle £2 p j x = An/N v \ x . Examples of such pixelisation schemes 
are HEALPix 1 iGorskietal. 1999) and IGLOO (Crittenden 2000). 
The distribution of pixel centres across the sphere is unimportant. 

Our simulation method begins by drawing each pixel value 
s p = S(x p ) — S(Qp,§p) independently from the same one- 
dimensional non-Gaussian PDF fs(s). The precise PDF used is 
unimportant, but for the purposes of illustration, we adopt here a 
PDF derived from the Hilbert sp ace of a linear harmonic oscilla- 
tor, as developed by Rocha et al. 1 2001 ). This class of PDF is sum- 
marised in Appendix A. In particular, we assume the PDF illus- 
trated in Fig.Q which is chosen for convenience to have a mean of 
zero. 

The resulting map S(x) will, by construction, be statistically 
isotropic to within the errors introduced by the pixelisation scheme. 
In fact, it is a realisation of isotropic non-Gaussian white noise, 
with a variance given by the second (central) moment jj 2 of the 
PDF fs(s). We will assume throughout that the mean /ji of the 
generating non-Gaussian PDF is zero, so that moments and cen- 
tral moments coincide. According to the 'cumulant expansion the- 
orem' iMa 1985) the cumulants, or connected moments, of the mul- 
tivariate distribution of pixel values are related to the logarithm of 
the moment-generating function, 



M(k) = (e*'*}, 



(1) 



where 5 = (si,.. .,s^ .J is the vector of pixel values and k = 
(/q , . . . , /c/y ix ) • The connected moments are given in terms of the 
derivatives of lnM(fc) by 



http://www.eso.org/science/healpix/ 
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Non-Gaussian Simulation 



Figure 1. The non-Gaussian PDF used in the simulations (solid line), which 
has the form given in equation 1A3I with (X3 = 0.2 and Oq = 1, and a set of 
samples drawn from the PDF (histogram). 
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Given that the pixel values s p are independent random variables we 
have lnM(k) = £ p lnMs(&p), where 



lnM s (k) = ln(e M ')=ln|~/ 5 (,) 



e isk ds 






(3) 



is the logarithm of the moment-generating function (the cumulant- 
generating function) for the non-Gaussian PDF fs(s) whose cumu- 
lants are the K„. Substituting in equation J2j> we find 

d" I 

( s Pi---Spn) c = Spv--P»(-i) n -^teM s {k)\ k=o = 8 Pl ... Pn K n . (4) 

The symbol 8 pl ... Pll equals one if all the n pixels are the same and 
vanishes otherwise. If required, the pixel correlations can always 
be expanded in their connected parts I Ma 1985). In particular, since 
(s p ) = for each pixel, one finds that, for example, 
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(Sp 1 Sp 2 Sp 3 Sp i jc + (S Pl Sp 2 )c{Sp 3 Sp i )c 
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The next step in the simulation procedure is to transform the 
map S(x) into spherical-harmonic space (using, for example, the 
map2alm routine from the HEALPix package) to obtain the coeffi- 
cients 



./47I 



Y ( * m (&,$)s(e, 40 do. 



(5) 



a lm — 2^i *fm\ x PJ s P*^pi]<- '' 
p=\ 

Using equation ^4) and the (approximate) orthogonality of the (pix- 
elised) spherical harmonics, we quickly find that the second-order 
correlator of the harmonic coefficients is given by 

{a( m a} m ,) =M2&pix&U'&mm', (6) 

where /J2 = K2 since fg(s) has vanishing mean. In order to obtain 
a final non-Gaussian map with a particular prescribed ensemble- 
average power spectrum, Q, one then rescales the harmonic coef- 
ficients to obtain 



a in 



Q 



/j 2 Q., 



(7) 




Figure 2. A realisation of a non-Gaussian all-sky CMB map with a pre- 
scribed ensemble-average power spectrum Q, obtained using the non- 
Gaussian PDF plotted in Fig.Qwith (X3 = 0.2, HEALPix resolution param- 
eter A/^ =512 (WMAP resolution) and C ax = 1500. 



Figure 3. A realisation of a Gaussian all-sky CMB map drawn from the 
same ensemble-average power spectrum, Q, as the non-Gaussian map 
shown in Fig. [2] 



such that (d( ! m S*f, m i) = Qo^f8 mm <. We note that the effect of a 
spatially-invariant, circularly-symmetric observing beam on the fi- 
nal map is trivially included by letting Q — ► QB?, where 6/ are the 
coefficients of the beam in a Legendre expansion. Finally, the har- 
monic coefficients d( m are inverse spherical harmonic transformed 
(using, for example, the alm2map routine from the HEALPix pack- 
age) to obtain the final map 



' = T{ x p) = lli a emYlm(Xp 



(8) 



Lin 



■pix 



where the double summation extends from £ = to £ max and m = 
—£ to I. The equivalent flat-sky approximation for small patches is 
discussed in Appendix Icl 

In Fig.|2|we plot a realisation of a non-Gaussian all-sky CMB 
map generated as described above, using the non-Gaussian PDF 
plotted in Fig. Q with a prescribed ensemble-average power spec- 
trum, Q. The map was produced using the HEALPix pixelisation 
scheme with the N s \i e parameter set to 512, which corresponds 
to Apix = 12A?? de ~ 3 x 10 6 equal-area pixels. For comparison, in 
Fig. [3] we plot a realisation of a Gaussian CMB map with the same 
power spectrum Q, using the same pixelisation. 

The source code to simulate the non-Gaussian CMB maps for 
both the full sky and for a small patch of the sky are available at the 
NGSIMS webpage 2 . 

As a guide to help the reader reproduce our method, we give 
here a summary of the logical steps to be taken to generate these 
non-Gaussian maps (see also documentation at the NGSIMS web- 
page): 

(i) Draw independent identically-distributed pixel values from 
a non-Gaussian PDF to create a statistically-isotropic map of non- 
Gaussian white noise; 

(ii) Transform the map to harmonic space (with e.g. a fast spher- 
ical transform); 



http://www.mrao.cam.ac.uk/~graca/NGsims/ 
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Power spectrum 




Figure 4. The power spectrum of the map shown in Fig.|2](points) as com- 
pared with the input ensemble-average spectrum Q (solid line). 



(iii) Scale the harmonic coefficients to enforce the desired power 
spectrum; 

(iv) Transform back to real space to obtain a non-Gaussian map 
which has the desired 2-point correlation function. 

We have implemented this method using two classes of non- 
Gaussian PDF (see documentation in NGsims webpage), but for 
the purposes of illustration, we adopt here PDF (i): 

(i) A general PDF based on the energy eigenstates of a linear 
harmonic oscillator, which takes the form of a Gaussian multiplied 
by a series of Hermite polynomials (see Appendix A). In princi- 
ple this expansion can be used to generate any PDF, although if 
the expansion is truncated then the available values of the relative 
skewness are constrained. 

(ii) The pixel values are drawn from a Gaussian distribution and 
then raised to some (even) integer power. This method is very fast 
to implement and easily generates distributions with large skewness 
and so is useful for checking statistical tools. 

As we shall see in the next section, the simulation method 
described above enables one to generate non-Gaussian maps for 
which many of the statistical properties can be calculated analyti- 
cally. The range of possible correlators and polyspectra are, how- 
ever, rather restricted, with the scale dependence of the latter con- 
trolled solely by the angular power spectrum (see Section l3~2t . It 
is, however, straightforward to extend our basic method to allow 
the simulation of non-Gaussian maps with a much wider range of 
statistical properties, as will be discussed in Sectional 



3 STATISTICAL PROPERTIES OF THE MAP 

It is clear that, by construction, the non-Gaussian map plotted in 
Fig. [2] is statistically isotropic, with an ensemble-average power 
spectrum, Q, corresponding to a generic inflationary cosmological 
model. As a check, the power spectrum of the map was calculated 
and found to agree with the input spectrum, as shown in Fig. [4] 
Owing to the simple manner in which the simulated non-Gaussian 
map is produced, many more of its statistical properties may be 
calculated analytically, such as the correlators of the pixel values, 
the bispectrum and higher-order polyspectra, and the marginalised 
probability distribution of the map. In this section, we again con- 
centrate on the all-sky case; the corresponding discussion for the 
flat-sky approximation is given in Appendix Icl 



To facilitate the calculation of the statistical properties of the 
map, it is convenient first to express the temperature t p in each pixel 
in terms of the initial pixel values {s p } drawn from the original 
non-Gaussian PDF fs(s). Using equations Q an d JSJ, we begin by 
writing 

: 2^ V Q«im^m(*p) 



V^^j 



■pi* t,m 



^ E v E n/qE^^M, 



(9) 



where in the second line we have substituted for ag m from equation 
J5J. Using the spherical harmonic addition formula 



2£+ 1 



4n 



(10) 



where P((z) is a Legendre polynomial, we may thus write equation 
{9} as 



V^E^VV' 
p' 

where the elements of the weight matrix are given by 



>w -^E^v^fe •*,')• 



(id 



(12) 



It is also useful to express the spherical harmonic coefficients, 
a£ m , of the map as linear superposition of the original pixel values 
{s p }. From equations J5J and 0, we immediately obtain 



P 



(13) 



where the elements of this second weight matrix are simply scaled 
spherical harmonics evaluated at the pth pixel and are given by 



W em , p = J^£YL(x 



K 



e.m\*p> 



(14) 



Equations il li — I14i provide the basic expressions from which 
the statistical properties of the non-Gaussian map may be obtained. 



3.1 Correlators of the pixel values 

Using equations ^4) and i\\\ . a general expression for the con- 
nected 72-point correlator of the pixel values in the non-Gaussian 
map is given by 

{'pi ■■■ t P,:)c = E W P< <H " ' W Pn9n Wl ' ' ' S 1n)c 

qi-q n 

= **Y,W M — W M , (15) 

q 

where K„ is the nth cumulant of fsis). Inserting the expression I12t 
into this result, we obtain 

K„ / O. ix \ "I 2 
(^.•••^„)c = 7T- , E a 4 •••4 / 4 •••4 fei. ••->*/>,, ),(16) 

"pix V M2 J iv .. (n 
where we have defined the quantities 

(2«i + !)••• (24 + 1) 
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V c ti-Ce., (17) 
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1 
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Taking the continuum limit of the sum over q, we see that an ana- 
lytic expression for the n-point correlator may be obtained by eval- 
uating 



h l -( / „( x pi'---' x p»)~ / p ti(*pi ■ x )--- p i„( x P„- x ) 

,/4it 



■x)d£l. 



(19) 



Clearly, h i ...i n is invariant under rigid rotations of its vector ar- 
guments, and this property is inherited by the n-point correla- 
tor as required by statistical isotropy. The integral may be re- 
lated to_the_w ;; p_olar harmonics with zero total angular momen- 
tum IVarshalovich et a l. 1988); this reduction for the first few val- 
ues of n is described further in Appendix 151 

Using the results derived in AppendixlBl one recovers, for ex- 
ample, the well-known result for the 2-point correlator 



Wpi) = L^- c e p (( x Pi ■Xpz)' 

and an explicit form for the 3-point correlator given by 

Vp\tpitpi)c 



(20) 
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in which we have defined Df = (2£ + 1)Q. 



3.2 Bispectrum and higher-order polyspectra 

To compute the polyspectra of the non-Gaussian map we form the 
nth connected correlator of its multipoles using equations ^4} and 
HI: 

{<3fimi ' • ' "£„m„ }c = 2- W dmi ,pi ' ' ' W e„m, u p„ ( s p, ' ' ' s p„ }c 

Pl-Pn 

= K„£w fimi , p ---ty />w ,. (2i) 

p 

Inserting the expression I14i into this result, we obtain 

K„ / SI i x \ "I 2 

{d( m l ---d( m ) c = — — \/Q, • • • Q„ hum ,.,4m, i (22) 

£i pix V K2 / 

where 

him,,...l„m„ = Z^^l! l m l ( x p)"'^t„m„( x p)^ i pi 



pix 



J 4k 



(23) 



Integrals of this form may be reduced to products of 3 j symbols as 
discussed in AppendixlBl 

For the case n = 2, one immediately recovers (a( m a% m ,} = 
C(8(p8 mm i. For n = 3, we find that 



£\ £2 £3 
(ae im a limi a hm ) = [ ^ ^ ) B ii4'« 



(24) 



where the bispectrum coefficients of the simulated non-Gaussian 
map are given by 
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(25) 



It is convenient also to introduce the 'normalised' (by the power 
spectrum) reduced bispectrum 8^4^ by the relation 
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(26) 



Thus, for our simulated map, b^.i, has a constant value given by 



b, 



K3 



tl*2«3 



Q, 



pix 



Q. \ 3 / 2 
i"2 



(27) 



We see that the amplitude of the bispectrum is determined by 
the variance /J2 an d skewness K3 of the original non-Gaussian PDF 
fs{s). If fs(s) were Gaussian, for example, the bispectrum would 
clearly vanish. Even if fs(s) has a non-zero skewness, however, 
we note from equation 1221 that, in the limit that the number of 
pixels tends to infinity, all polyspectra with n ^ 3 vanish. This is 
a consequence of each pixel being the weighted sum of indepen- 
dent identically-distributed variates, so that in the limit of an infinite 
number of pixels the processed map T(x) tends to Gaussian by the 
central-limit theorem. Fortunately, in practice, we can bypass this 
property as /V p i x increases by scaling the cumulants of the initial 
distribution fs(s). For example, to get a bispectrum independent of 
Afpjx, we must scale K3 such that 



v pix 

fin 



Q. \ 3 / 2 

"pix 



= constant. (28) 

'pix \ A<2 

As mentioned in AppendixlAl for the particular non-Gaussian PDF 
used here, generating higher values of the relative skewness re- 
quires one to have a larger range of the generalized cumulants pa- 
rameters a„ non zero. 

As a check on our calculations, the normalised reduced bis- 
pectra of an ensemble of 300000 non-Gaussian maps was calcu- 
lated using an estimator defined in Sperael & Goldberg 1 1999) as 
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(29) 



where, ei(x) = \/ jfrjY.m a lmYl'm( x )- The resulting mean values of 
bin are plotted in Fig.|5] together with the associated uncertainties. 
Note that it is important to divide the measured bispectrum by the 
ensemble-averaged values of the power spectrum; dividing by the 
measured values for each simulation produces a biased estimator. 
The predicted value of P44& was calculated using 1271 and is plot- 
ted as the solid line in the figure. We see that the measured and 
predicted values are fully consistent. 

Returning to equation \12\ . if one considers n = A and follows 
the notation of Hu 1 2001), one finds that the connected part of the 
trispectrum of the non-Gaussian simulation is given by 



P lt( L ) = 



f'pixK4 
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yjD tl D t2 D l3 D u 

' l\ h L 




2L+1 

An 



£ 3 £4 L 




(30) 



Analytic expressions for higher-order polyspectra may be obtained 
in an analogous manner. 

3.3 Cumulants and marginal distribution 

Finally, we consider the marginal PDF fr(t) of the processed map. 
This is mostly easily carried out by considering the cumulants of 
fr(t), and relating them to those of the original non-Gaussian PDF 

fs(s)- 
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Full-sky non-Gaussian simulation 



o 



0.035 

0.03 

0.025 

5- 0.02 

f 0.015 
j 

0.01 
0.005 















Measured values i — i — > 












-i-l--i-L-i-j-I-i-^-T-i-T-i-T-i---i-i-" 


- 






\i] k yi i jppnii *_ 




Figure 5. Non-zero diagonal components of the bispectrum estimated from 
300000 non-Gaussian simulations with the same parameters as those used 
in Fig.|2| The mean over the simulations and its standard error are plotted. 
The dashed line is the theoretical ensemble-average value. Note that the 
diagonal bispectrum necessarily vanishes for odd I. by parity. 



The cumulants of the marginal PDF are equal to the connected 
parts of the correlators of the pixel temperatures with all pixels in 
the correlator the same. From equation 1151 . we find that 
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which is independent of the choice of pixel p. Evaluating the 
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axis, where Y( m (z) = y/(2l + l)/47t5 m n, we find that the first few 
cumulants (beyond k'j = 0) are 
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where we have written the results in such as way that they are true 
generally, for any statistically-isotropic map. 

In principle, knowledge of the complete set of cumulants K r n 
may be used to obtain an explicit expression for the marginal PDF 
fr{t)- This could be carried out, for example, by first obtaining its 
moment-generating function 



M r (fc) = (e lfa )=exp(£^K'„ 

\n=l "• 



(32) 



and then performing an inverse Fourier transform to yield frit). 
Alternatively, for a weakly non-Gaussian distribution, one can em- 
ploy the Edgeworth expansion. In this approach, the PDF is ex- 
pressed as an asymptotic expansion around a Gaussian with mean 



Figure 6. Histogram of pixel temperatures for the final processed non- 
Gaussian map for A/ SK / C = 64 and £ max = 128, created using an initial PDF 
with parameters 013 = 0.27 and o"o = 1 . Also shown are the Edgeworth ex- 
pansion for fr(t) (dashed line), given by equation 1331 but truncated at 
n = 4, and a Gaussian with the same variance (solid line). 
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where t — t/yTxs. When kJ,/o" is small for n larger than some in- 
teger one can use a finite number of cumulants as an acceptable 
approximation to the distribution. As pointed out by Roc ha et alJ 
(2001), however, this might no longer be a PDF and, in particular, 
might deviate from the original distribution in the tails. In Fig|6| 
we plot the histogram of the pixel temperatures for a non-Gaussian 
map with A^ide = 64 and £ max =128 (somewhat lower resolution 
than that plotted in Fig.|2j, generated from an initial PDF with pa- 
rameters 0C3 = 0.27 and Oq = 1. Overplotted is the Edgeworth ex- 
pansion of fr(t), equation I33t . truncated at n = 4. The cumulants 
of frit) can be computed efficiently from the first expression in 
equation J3D; we find k\ = 0, K^ w 13809 pK 2 , k£ w 177234,uK 3 
and K^ ~ 6307963 ,uK 4 . The Edgeworth expansion agrees with the 
simulation results better than a Gaussian with the same variance. 



4 EXTENDED SIMULATION METHOD 

We see from the previous section that our basic simulation method 
generates maps with a rather restricted range of possible correlators 
and polyspectra, with the scale dependence of the latter controlled 
solely by the angular power spectrum. It is straightforward, how- 
ever, to extend our basic method to allow the simulation of non- 
Gaussian maps with a much wider range of statistical properties. In 
particular, there is no fundamental requirement for the method to 
be restricted to the same set of cumulants over the whole range of 
scales. 

The procedure is as follows. One divides the range of mul- 
tipoles, I, into non-overlapping bins. A given bin B will be a set 
B ~ [^Bmini^max]- ^" 0T eacn bin B we simulate a map tgp, with non- 
zero Q only for £ £ B and a non-Gaussian PDF that can differ be- 
tween bins. The final map is a superposition of these band maps, 
i.e. with pixel values 
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(34) 



Simulation of non-Gaussian CMB maps 



Non-Gaussian Simulation usine 3 bins 




O 



0.03 - 
0.025 
0.02 

4. 0.015 -_ 

a" 
v 

0.01 - 
0.005 

- 



Full-sky non-Gaussian simulation using 3 bins 



Measured values 
Theoretical values 



Figure 7. A realisation of a non-Gaussian all-sky map simulation using 
three bins: B { = [0,500], a 3 = 0.1; B 2 = [501, 1000], a 3 = 0.2; and B 3 = 
[1001, 1500], a 3 =0.25. 



Since each map tg p is individually statistically isotropic, then so 
too is their sum. Moreover, from equation <26> . we see that the 
bispectrum for each band map can only be non-zero within the cor- 
responding bin (this is also true for the connected parts of higher- 
order poly spectra). The «th cumulant of the final map is also simply 
the sum of the nth cumulants of the individual band maps. With this 
method we are thus able to generate maps with more general statis- 
tical properties. 

For example, by choosing appropriate values of K3 for the non- 
Gaussian PDF used to simulate each band map, one can arrange for 
the summed map 1341 to have a given power spectrum Q and an 
arbitrary prescribed constant value of the reduced normalised bis- 
pectrum £%^ 3 in each bin. As an illustration, in Fig. we plot 
a non-Gaussian map generated using three bins. In each bin, the 
non-Gaussian PDF used was of the form given in equation <A3> 
with CTo = 1. However, the values of 0C3 used in each bin were 
(X3 = 0. 1 , 0.2 and 0.25 respectively. As a check on our calculations, 
the normalised reduced bispectrum of an ensemble of 300000 such 
non-Gaussian maps was calculated using the estimator 1291 . The 
resulting mean values of hue, for individual values of £, are plotted 
in Fig. [5] together with the associated uncertainties. The predicted 
value of bin m eacn of the three broad bins was calculated using 
equation 1271 and are plotted as the dashed lines in the figure. We 
see that once again the measured and predicted values are fully 
consistent. 

How general is our extended method? We can in principle ef- 
ficiently generate maps with arbitrary diagonal bispectra, i.e. with 
any given Q and Bfn. This is of some importance because cur- 
rently known primordial theories of non-Gaussianity lead to more 
general combinations of angular power spectra and bispectra than 
can be created from our method with a single univariate PDF fs(s). 
It is not possible, however, to generate specific models of primor- 
dial non-Gaussianity exactly with our extended method. To do that 
one would have to be able to choose arbitrary values for the angu- 
lar spectrum and for all components of the bispectrum (and higher 
polyspectra). This would involve going beyond the simple one- 
point PDF methods advocated here. 



5 CONCLUSIONS 

We presented a simple, fast method for simulating statistically- 
isotropic non-Gaussian CMB maps with a given power spectrum 



Figure 8. Non-zero diagonal components of the bispectrum estimated from 
300000 non-Gaussian simulations using three bins: B\ = [0,21], W3 = 0.1; 
B 2 = [22,43], a 3 = 0.2; and B 3 = [44,56], a 3 = 0.25. The mean over the 
simulations and its standard error are plotted. The dashed lines shows the 
theoretical ensemble-average value in each bin. 



and analytically calculable bispectrum. We showed that our tech- 
nique allows one to describe the statistical properties of the map 
by computing analytically the nth-order polyspectra, and the nth- 
order correlators of the pixel values. We showed that these can 
be expressed in terms of the rc-polar harmonics with zero total an- 
gular momentum, and we describe this reduction for the first few 
values of n. We also recovered analytically the one-dimensional 
marginalised distribution function in terms of its cumulants. The 
univariate non-Gaussian distribution, from which the pixel values 
are drawn independently in the first stage of the simulation process, 
fully determines the statistical properties of the final map. Here we 
used a non-Gaussian distribution derived from the wavefunctions 
of the harmonic oscillator. Simulations of both the full sky and a 
small patch of the sky were generated and corresponding statistical 
analysis performed. As a check on our calculations we computed 
both the power spectrum and bispectrum of the simulated maps and 
found them to be fully consistent. The simulation method described 
here clearly enables one to generate maps with well-defined corre- 
lators and polyspectra. We extended the method to encompass dif- 
ferent set of cumulants over the whole range of scales, generating 
maps with arbitrary power spectra and diagonal bispectra for differ- 
ent scales. It is not possible, however, to generate specific models 
of primordial non-Gaussianity exactly with our extended method, 
since these require off-diagonal bispectrum coefficients to be spec- 
ified arbitrarily. This would involve going beyond the simple one- 
point PDF methods advocated here. 

The source code to simulate the non-Gaussian CMB maps for 
both the full sky and for a small patch of the sky are available at the 
NGSIMS webpage 3 . 

A pertinent question is what other statistical properties can be 
calculated analytically for the class of non-Gaussian maps we have 
investigated. Of particular interest are the phase associations be- 
tween different harmonic coefficients. In Matsubara 12003), a gen- 
eral relationship between phase correlations and the hierarchy of 
polyspectra in Fourier space is established. It is also stated that the 
phase correlations are related to the polyspectra through the non- 
uniform distribution of the phase sum 6^ + 0^ H + % with 

closed vectors k\ + &2 H +^A< = 0. We are currently investigat- 



http://www.mrao.cam.ac.uk/~graca/NGsims/ 
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ing the form of the distribution function of this phase sum in our 
maps. A study of the Minkowski functionals of our non-Gaussian 
maps is also underway. 



Vio R., Andeani P., Tenorio L., Wamsteker W., 2002, PASP, 114, 
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APPENDIX A: NON-GAUSSIAN PDFS BASED ON THE HARMONIC OSCILLATOR 

In this appendix, we summarise th e class o f probability distribution functions (PDFs) derived from the Hilbert space of a linear harmonic 
oscillator, which was developed bv lRocha e t al. 1 2001). The original non-Gaussian distribution, fs(s), used in the main text to produce the 
simulated non-Gaussian maps is an example of such a PDF. 

This general PDF is based on the coordinate-space wavefunctions of the energy eigenstates of a linear harmonic oscillator, and takes 
the form of a Gaussian multiplied by the square of a (possibly finite) series of Hermite polynomials whose coefficients rx„ are used as 
non-Gaussian qualifiers. In particular, if x is a general random variable, the most general PDF has the form 



P (x) = \m 2 



-*7(2og) 



2^a„C n H„l —= — 
„ \ V2oo 



(Al) 



where H n {x) are the Hermite polynomials, and the quantity Oy is the variance associated with the (Gaussian) probability distribution for the 
ground state \\p$\ ■ The constants C n are fixed by normalising the individual states. The only constraint upon the amplitudes rx„ is 

£|ce„| 2 = l. (A2) 



This is a simple algebraic expression which can be eliminated explicitly by writing oCq = Jl — £" |oc„| 2 . Thus the coefficients oc„ can be 
independently set to zero without mathematical inconsistency iRoch a et alJl200ll) . Moreover, these coefficients can be written as series of 
cumulants (Contaldi. Bean & Maaueiio 1999) and should indeed be regarded as non-perturbative generalisations of cumulants. 

For the simulations in the main text, we use the non-Gaussian PDF for which all a„ are set to zero, except for the real part of CC3 (and 
consequently do). The reason for this choice is that this quantity reduces to the skewness in the perturbative regime. The imaginary part of 
CX3 is only meaningful in the non-perturbative regime (and can be set to zero independently without inconsistency). Hence we consider a PDF 
of the form 



p(x) 



B -*?/(2o5) r a 3 2 



2no<\ 



a + -/=#3 



V48 V\/2o 



(A3) 



with otp = \/ 1 — a^. It is straightforward to show that the first, second and third moments of our PDF are related to 0:3 and Cq by 
JContaldi & Ma P ueiiol200lh 

A/2 = 0^(1+60^ 



W = (2cly ^J3[al(l-ai)}. (A4) 

The PDF therefore has zero mean and a fixed variance and skewness. In the simulations discussed in the main text, we choose 0:3 = 0.2 and 
o"o = 1. This resulting PDF is plotted in Fig.Q 

We note that the space of possible PDFs is constrained as a result of restricting the set of coefficients a„ to two non-zero values. This 

3/2 
implies that we cannot generate distributions with arbitrarily large relative skewness. Indeed, /J3/1U2 is bounded above by 0.74, and takes 

this maximum value for (X3 = (7 — \/43)/6 = 0.27 . However, in general our method can generate highe r values of the relative skewness 

(since it can generate any distribution) but for that purpose one needs more non-zero coefficients oc„ (Contaldi & Magueijo 2001). 

APPENDIX B: SOME USEFUL INTEGRALS 

We give here useful results concerning integrals involving products of Legendre polynomials: 

k v .. ln {x x ,...,x n )= J P tl (x x -x)---P ln {x n -x)&Q.. (Bl) 

Using the addition theorem (equation llOt . this reduces to evaluating integrals of products of spherical harmonics, 

Jim, e„ m „ = / Ye l m l ( x )--- Y e„m„(x)dQ. (B2) 

JA% 

since 

471 471 

k w ..l n {x u ...,x n ) = ——■■■—— £ Yl (xi)---Yl mn (x n )Jt i m u ...Am n - (B3) 

First, consider the case n = 2. Using the orfhonormality of the spherical harmonics, and the relation Yi L(x) = (—l) m Y(_ m (x), to evaluate 
J( mi ( 2 „, 2 = (— l) mi 8^ 2 Smi -'"2' anc ' tnen a PPly m g me addition theorem we find the well-known result 

47t 

hii 2 (*uX2) = 2g rj -fyfa -x^uh- ( B4 > 

For integrals involving product s of three or more s pherical har mon ics, the general strategy is to combine pairs of harmonics using the 



^or integrals involving products ot three or more spherical harmon: 
Clebsch-Gordan series (e.g. lVarshalovich et all 1 988: Edmonds 1974) 
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Y, MY, M-V, / (^l + l)(^2 + l)(2^+l) / £ l £ 2 l\(ll h t\y* (x) ms) 

^ mi (m 2m2 ^)-L\J ^ { m m m )\o o o ) Y ^ X >' (B5) 

until we have only a single pair left which can then be integrated trivially using orthonormality. We illustrate this procedure for the case 
of n = 3 and n = 4 since these are needed for the calculation of the bispectrum and trispectrum. For n = 3, from equation IB5i we find 
immediately that 



/ (2l 1 + l)(2l 2 + l)(2l3 + l) f h h h \( h h h i ,„,, 

Jl imi ,...l imi -\J ^ ^ )[ mi mi m ). (1 > ( » 

and therefore 



■ 2 / *X Ml <2 <3 A r- / <1 ^2 ^3 



/ -"'— Mit ^/ (2. 1+ i)(2£ 2+ i) ( 2.3 + i) U o oJ m ,LU ** m ;j 7 ^ (jci)F ^ ( ^- fe) - (B7) 

The final term in this equation (the summation over m\, m 2 and m^) ensures that hit z t^ ls invariant under rigid rotations of its vector 
arguments x \ , xj, and xt,. As expected , the summation can be expressed in terms of the tripolar spherical harmonics with zero total angular 
momentum l Varshal ovich et alll 988). 

Consider now the case n = 4. There is now some freedom in the choice of spherical harmonics to combine. If we couple Y( imi with Y# m2 
and y> 3 ,„ 3 with Y lim , we find 



T V( IV » / (2£i + l)(2£ 2 + l)(2£+l) / (2€ 3 + l)(2€ 4 + l)(2€+l) 

' <i ^2 M ( h U i \( h h l \f h k i \ (B9) 

m\ m 2 tn J \ mj, 1114 —m J \ OyyO J 

The expression on the right is not manifestly symmetric with respect to interchange of e.g. (l\m\) and (£31112) since the latter in- 
volves a different coupling scheme. However, the symmetry is easily verified by switching between the two schemes with the 6j coeffi- 
cients 1 Varshalo vich et all 19 88V Finally, we find that 



,a n3/2 /(2£j + l)---(2f 4 + l) v- 
I lv .. h {x h ...,X4) = (4K) il2 \j± '-^ L Y, 



(ai+i)-(a 4 +i) r 2i+i/<i £ 2 £ \ (£ 3 £ 4 £ 
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It is straightforward to verify that the last term on the right (the summation overra, mi,... ,1114) is invariant under rigid rotations ofjci,. .. ,xa,, 

APPENDIX C: FLAT-SKY APPROXIMATION 

For analysis over a small patch of the sky we can use the flat-sky approximation and replace spherical transforms by Fourier transforms. 
Our starting point is again a pixelised map of non-Gaussian white noise, with each pixel value drawn from the non-Gaussian PDF fs(s). 
Approximating the Fourier transform a(l) by a discrete Fourier transform we have 

a{l) = j^S{x)s-** ~ Sej^e-**', (CI) 

where Q. p [ x is the pixel area. We evaluate a(£) on a regular grid in Fourier space with a Fast Fourier Transform. For a square patch of sky 
with A'pjx pixels, the cell size in Fourier space is (2k) / (N p [ x Q, p [ x ) . The second-order correlator of the discrete a(£) evaluates to 

"pix 



(a(£)a*(£!)) = l-^\ N&i$ u , (C2) 

where /j 2 is the variance of the zero-mean f$(s). In the continuum limit, equation IC2i becomes 

(a(e)a*(e , ))=Op xM2 s(e-e'), (C3) 



where we have used 

\2 



8 = J_£ Jft-i)*, _> _JL_ [# x j<!-*)* = J^8(£-t). (C4) 

iv pix p ^pix^pix ■-' ^ v pix"pix 



We scale the a(£) defining a(£) = ^/Q/(£2 p i XJ u 2 )a(^) (as for the full-sky case described in the main text) such that the a(£) have the required 
power spectrum: 

{d(e)a*{e'))=c ( 8(£-e'). (C5) 

Finally, we inverse Fourier transform to obtain our non-Gaussian map, T(x), with the prescribed two-point statistics: 
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T (x p ) = -J^—y^y lx " « f —a(ey lx " ■ (cq 



271 v -,„s if.r f d £ . 

v pix ii pix 

As in the full-sky case, we can express the pixel values, f p , in the final map as linear combinations of those in the original map, s p : 



p' 

where in the continuum approximation 



"s* - s^/y^'""-'" 



— J^VQJoWxp-Xp>\), 



with 7o(z) the Bessel function of order zero. Using the asymptotic result Jo(£Q) ~ P/>(cosQ), it is straightforward to see that W pp i obtained 
here in the flat-sky limit is equivalent to the full-sky expression (equation 1 12> . Forming the connected rc-point function, as in equation I15i . 
we find 

For example, for n = 2 we have 

f d 2 £ \f(rr\ t Idl 

(t P M)c = J T^iQe * ( *" n) « / ^■QJbW*?, -**l)- (cio) 

In the asymptotic limit this result reduces to the full-sky expression (equation l20t . 

Finally, we consider the polyspectra of the processed non-Gaussian maps. Expressing the Fourier transform of the final map, a(£), in 
terms of the original map, i.e. 



M-iJ^Z'p*-**'' < C11 ) 
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^L^fVy (27t) 2 -VQ 1 -Q,8(< 1 +•••+/„), (C12) 

where in the last line we have taken the continuum approximation. This form for the correlator is clearly consistent with rotational, transla- 
tional and parity invariance. 

As in the full-sky case, we have produced simulated non-Gaussian maps and calculated their power spectra and bispectra; the latter were 
estimated using the code described in lSmith e t al. 12004). We find that the flat-sky simulations behave as expected with no discernible bias. 

This paper has been typeset from a TpX/ KTpX file prepared by the author. 



